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Toward Bernal Random Loose Packing through freeze-thaw cycling 
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We study the effect of freeze-thaw cycling on the packing fraction of equal spheres immersed in 
water. The water located between the grains experiences a dilatation during freezing and a con¬ 
traction during melting. After several cycles, the packing fraction converges to a particular value 
rjoo = 0.595 independently of its initial value r/o. This behavior is well reproduced by numerical 
simulations. Moreover, the numerical results allow to analyze the packing structural configuration. 
With a Voronoi partition analysis, we show that the piles are fully random during the whole process 
and are characterized by two parameters: the average Voronoi volume jj, v (related to the packing 
fraction rj) and the standard deviation a v of Voronoi volumes. The freeze-thaw driving modify 
the volume standard deviation a v to converge to a particular disordered state with a packing frac¬ 
tion corresponding to the Random Loose Packing fraction t/brlp obtained by Bernal during his 
pioneering experimental work. Therefore, freeze-thaw cycling is found to be a soft and spatially 
homogeneous driving method for disordered granular materials. 

PACS numbers: 81.05.Rm,45.70.Cc 


How a large number of solid objects can fill a vol¬ 
ume is one of the most puzzling problems in mathemat¬ 
ics, science and engineering. This question concerns a 
broad range of systems: granular media, colloids, struc¬ 
tures of living cells, amorphous solids, ... The relevant 
parameter that characterizes a pile of particles is the di¬ 
mensionless packing fraction 77, defined as the volume 
of all particles divided by the apparent volume of the 
assembly. This packing fraction has a maximum value 
rjfcc = tt/ 3%/2 ~ 0.74 for identical spheres, corresponding 
to the face-centered cubic (fee) lattice. The existence of 
this maximum was stated by Kepler in 1611 and demon¬ 
strated recently Tj. The first systematic study of random 
equal spheres packings has been performed by Bernal et 
al [2|. The random packing of spheres was considered as 
a useful model for ideally simple liquid. The obtained 
range of random packing fraction was found to lie be¬ 
tween two well-defined limits: the Bernal Random Loose 
Packing fraction tjbrlp = 0.60 and the Bernal Random 
Close Packing fraction tjbrcp = 0.63. The value of the 
packing fraction depends strongly on the history of the 
pile, i.e. on the way particles are sequentially placed in 
the assembly. More recently, the range of packing faction 
for random packings has been extended. By using flu¬ 
idized bed techniques [3] or in a fluid providing a strong 
buoyancy [1J , one can decrease the random loose packing 
fraction to t]rrp ~ 0.55. The present accepted value for 
the random close packing limit is around rjRCp — 0.64 
[5] . The exact value of these limits and the link with the 
packing structure is still a matter of intense debate [SHE]. 
The situation is even more complex with non-spherical 
grains mm- 

Different driving methods are leading to the modifica¬ 
tion of the packing fraction: tapping HM, cyclic shear 
m , flow-induced fluidization BE •• . The existence of 
cohesive forces is also known to modify the packing frac¬ 
tion [16]. Recently, the effect of thermal cycling on pack¬ 
ings has been studied experimentally nans] and numeri¬ 
cally HE Despite the low magnitude of the induced ther¬ 


mal expansion, a significative densification of the packing 
has been observed. 

Beyond the classical problem of packings, many phe¬ 
nomenons observed in nature are related to the particular 
behavior of wet granular materials submitted to temper¬ 
ature cycling: ice-lens formation in soil leading to frost 
heaving mm, landslides, structures formation in per¬ 
mafrost, stone heave and possibly some geological for¬ 
mations observed on Mars. Indeed, both granular ma¬ 
terials and water have remarkable properties. The wa¬ 
ter in the liquid phase has a density maximum around 
4°C and the solid phase density is lower than the liq¬ 
uid phase density. Moreover, comparing the other ma¬ 
terials, water has a very high specific heat capacity, a 
high heat of vaporization, a high latent heat and a high 
surface tension. On the other hand, granular materials 
are out of equilibrium systems witch exhibits remark¬ 
able behaviors: phase segregation [221 . jamming transi¬ 
tion [23], extremely slow compaction dynamics |1H413) . 
glassy dynamics j26j, ■ • • Therefore, the association of wa¬ 
ter and granular material leads ineluctably to rich behav¬ 
iors [2i[ [25], 

In the present letter, we show how the packing fraction 
of a spheres pile totally immersed in water is modified by 
successive freeze-thaw transitions. Moreover, the evolu¬ 
tion of the packing structure during freeze-thaw transi¬ 
tions have been analyzed with Voronoi partition method 
applied on packings obtained by numerical simulations. 

A sketch of the experimental set-up is presented in 
Figure [l] A borosilicate glass tube of internal diame¬ 
ter D tu b e = 18.4 mm is filled with water. Afterward, a 
granular material made of glass spheres with a diame¬ 
ter Dg ra i n = 0.5 mm is poured gently inside the tube. 
This method leads to a loose initial packing fraction. To 
obtain a higher value of the initial packing fraction, the 
tube is placed inside an ultrasonic cleaning tank during 
one minute. The vibrations induce a compaction of the 
granular bed. To perform temperature cycling, the tube 
is placed inside an isotherm box. The temperature is 
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FIG. 1: (Color online) Sketch of the experimental setup. A 
high CCD camera records the top of the packing placed in 
a glass tube. The pile is back-illuminated by a homogenous 
lighting system. 
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decreased with a Peltier module. The temperature in¬ 
crease is produced by a resistive wire winded around the 
glass tube. Two thermocouples measured respectively 
the temperature of the air inside the box and the tem¬ 
perature inside the granular material. A CCD camera 
takes pictures of the top of the granular pile. The posi¬ 
tion of the water meniscus is also visible on the pictures. 
The height of both water/air and granular/water inter¬ 
faces are obtained by image treatment. Test runs with 
no thermal cycling were performed to ensure that vibra¬ 
tions do not affect the measurements. Moreover, thermal 
cycling have been performed on granular pile without wa¬ 
ter. In this dry case, we do not observe any significative 
variation of the packing fraction with our experimental 
set-up. 

The volume thermal expansion coefficient of the 
glass beads and of the tube are respectively /3 g i aS s = 
25.5 fiR- 1 and fiboro. glass = 9.9 fiK^ 1 . The water 
thermal expansion coefficient varies significantly with the 
temperature and ranges between —50 ^A ' _1 to 207 fiK -1 
when the temperature goes from 1°C to 20°C. A rough 
estimation of the different material expansion for a tem¬ 
perature ramp between 0°C and 20°C gives 0.02% for 
the tube, 0.05% for the grains and 0.15% for the water. 
Moreover, during the freezing transition, the water vol¬ 
ume increase is around 8 %. These estimations show that 
the presence of water in the system is expected to influ¬ 
ence deeply the thermal cycling compaction dynamics, in 
particular during freeze-thaw transitions. 

The pile is submitted to freeze-thaw cycling with tem¬ 
perature ranging from = —4 °C and T 2 = 24°C (see 
Figure [2] (bottom)). The pile temperature reaches neg¬ 
ative values before freezing. Therefore, a supercool¬ 
ing effect is observed. The time necessary to obtain 
a freezing of the system fluctuates from one cycle to 
the other. During the freezing transition, the temper- 


FIG. 2: (Color online) Top: picture of the top of the pile (a) at 
the beginning of the cycle, (b) just before the freezing, (c) just 
after the freezing, (d) at the end of the freezing procedure, (e) 
during the heating and (f) at the end of the cycle. Bottom: 
temporal evolution of the temperature inside the pile, of the 
water height and of the granular height during one freeze-thaw 
cycle (r/o = 0.573). 


ature inside the pile goes quickly to 0°C. During this fast 
freezing step, the position of the granular/liquid and liq¬ 
uid/air interfaces do not change significantly (see Figure 
[2] (Top)). However, during the temperature decrease af¬ 
ter this freezing transition, both interfaces are going up. 
Then, a dilatation of the whole system is observed. Some 
steps are observed during this increase. When the pile 
temperature reaches T\ = —4 °C, the Peltier module is 
switched off and the heating is powered on. As a conse¬ 
quence, the temperature inside the pile increases to the 
temperature T 2 = 24° C with a plateau at 0°C. During 
the heating, the interfaces are going down. Some steps 
are also observed during this heating process. 

The evolution of the packing fraction rj as a function of 
the cycle number n is presented in Figure [3] (top). Three 
experiments performed with different values of the initial 
packing fraction 770 are presented. The experimental data 
are well fitted by the single exponential law proposed by 
Mehta et al [27] v( n ) = Voo ~ A ye~ n ' T , where r/oo, Ay 
and r are respectively the asymptotic packing fraction, 
the range of packing fractions and a characteristic cycle 
number. For low initial packing fraction, the freeze-thaw 
cycling induces a densification of the pile ( 77 ^ = 0.592 ± 
0.001, A ?7 = 0.019 ± 0.001 and r = 3.5 ± 0.3). On the 
other hand, a decompaction is observed for higher initial 
packing fraction (7700 = 0.596 ± 0.002, A77 = —0.021 ± 
0.001 and r = 3.4 ± 0.7). In both case, the packing 
fraction is found to converge to the Bernal Random Loose 
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FIG. 3: (Color online) Evolution of the packing fraction p as 
a function of the freeze-thaw cycle number for experiments 
(top) and simulations (bottom). The squares and the circles 
are corresponding respectively to a loose and a close initial 
packing. The triangles are corresponding to an initial packing 
fraction close to the Bernal Loose random Packing fraction 
Vbrlp = 0.60. The curves are fits from exponential law (see 
text). 



FIG. 4: (Color online) Representation of the pile during one 
simulated freeze-thaw cycle. The grains are colored in yellow 
and the void spheres corresponding to the ice are colored in 
blue (color online). From the left to the right: initial pile with 
a low packing fraction, the void spheres are added, dilatation 
of the void spheres, final state with a higher packing fraction. 


Packing fraction tjbrlp = 0.60. Moreover, the packing 
fraction is not influenced by freeze-thaw cycling when 
the initial packing fraction is close to pbrlp- Therefore, 
this particular packing fraction should corresponds to a 
specific structural configuration of the packing. 

Although the compaction was expected for a loose ini¬ 
tial packing fraction, the decompaction process of a dense 
pile is surprising. Moreover, the convergence of the pack¬ 


ing fraction to the Bernal Random Loose Packing fraction 
iIbrlp = 0.60 is striking. In order to analyze the struc¬ 
tural evolution of the packing during freeze-thaw cycling, 
numerical simulations have been performed. The model 
is based on molecular dynamics with tangential spring in 
order to produce a static pile. To obtain details about 
the numerical method, see ref [3S]. A tube of diameter 
Dtube = 10mm is filled with N=1850 grains of diame¬ 
ter = 1mm. The grains have the density of glass 

(p = 2500 kg m~ 3 ). After the initialization method, the 
voids between the grains are filled with spheres having 
the density of water. The void positions are obtained 
by Voronoi' tessellation. Typically, 8000 void spheres are 
added. A freezing step consists in the dilatation of the 
void spheres, while a thaw step corresponds to a con¬ 
traction of the void spheres. The evolution of the pile 
during a simulated freeze-thaw cycle is presented in Fig¬ 
ure [4] When two void spheres are overlapping, a spring 
link is defined between them in order to avoid collapse 
during the dilatation process. As shown by Figure [3] 
(bottom), the compaction for loose initial packing frac¬ 
tion (r/ao = 0.594 ± 0.001, Ap = 0.019 ± 0.001 and r = 
9.7± 1.2) and the decompaction for higher initial packing 
fraction ( 77 ^ = 0.600 ± 0.001, A ij = —0.005 ± 0.001 and 
r = 5.0 ± 0.5) are well reproduced by the simulations. 
The characteristic cycle number obtained in simulations 
and in experiments are different because the size of the 
grains and of the container are different in experiments 
and in simulations. 

The packings obtained by numerical simulations have 
been analyzed with Voronoi partition. The probability 
density functions of the Voronoi cell volumes V for all 
the initial piles and for all the piles after 20 freeze-thaw 
cyclings are presented in Figure [5] The volume normal¬ 
ization (V — p, v )/cr v with the Voronoi volume average 
p v and the Voronoi volume standard deviation cr v in¬ 
duces a collapsing of the density functions. Then, the 
volume average p v and the volume standard deviation 
a v are the main parameters characterizing the system. 
The volume average p v is related to the packing fraction 
V = Vgrainl AG > where V gra i n is the volume of one grain. 
Moreover, the Gamma shape of the density functions is 
characteristic of a fully random system [29]. We have 
tested with numerical simulations that an other driving 
mechanism like vibrations induces a deviation from the 
Gamma shape due to the apparition of ordered domains. 
Therefore, freeze thaw cycling is found to be a soft driv¬ 
ing method for disordered granular materials. The driv¬ 
ing modify the volume standard deviation a v to converge 
to a particular disordered state with a packing fraction 
close to tjbrlp = 0.60. Contrary to driving processes 
based on mechanical agitation, we do not observe any 
convective motions, any grains spatial organization in the 
bulk and any grains ordering close to the container wall 
mm- Moreover, the driving method based on freeze- 
thaw cycling is spatially homogenous. 

A diagram with the main parameters obtained from 
the simulated packings (packing fraction p = V gra i n /p v 
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FIG. 5: (Color online) Probability Density Function (PDF) 
of the Vorono'i volumes V for the initial piles (open symbols) 
and after 20 freeze-thaw cycling (plain symbols). The vol¬ 
umes V are normalized with the volume average fj, v and the 
volume standard deviation a v . The convention for the symbol 
shapes are the same than in Figure [3] The data have been 
obtained from numerical simulations. The line corresponds 
to the Gamma function. 



FIG. 6: (Color online) Evolution of the Vorono'i volume stan¬ 
dard deviation a v as a function of the packing fraction r/ dur¬ 
ing 20 cycles for three different initial packing fractions. The 
data have been obtained from numerical simulations. The 
dashed curves are guide for the eyes inspired from the work 
of Aste et al. [6]. 


and standard deviation a v of the Vorono'i volumes) is pre¬ 
sented in Figure [ 6 j Aste et al. [ 6 ] have shown that a v 
decreases between the random loose and the random close 
packing limits. Moreover, they observed a small but siz¬ 
able local minimum around 77 = 0.60. This behavior was 
observed with experimental and numerical granular as¬ 
semblies created in different conditions. The freeze-thaw 
cycling allow to move in this diagram along the master 
curve describe by Aste et al.. We have checked numeri¬ 
cally that a stronger mechanical driving mechanism like 
vibrations induces a deviation from this master curve due 
to the apparition of ordered clusters in the packing and 
close to the wall. 

In summary, the evolution of the packing fraction of 
an assembly of spheres immersed in water and submitted 
to freeze-thaw cycling has been investigated. As already 
observed in the dry case, the packing fraction increases 
when starting from a loose configuration. However, when 
starting from a close configuration, the freeze-thaw cy¬ 
cling induces a decompaction of the pile. Independently 
of the initial packing fraction 770 , the packing fraction 
converges to a particular value 7)^ = 0.595. This fi¬ 
nal packing fraction is equivalent to the Random Close 
Packing fraction obtained by Bernal during his precur¬ 
sor experimental works performed in the sixties. Con¬ 
trary to driving processes based on mechanical agitation, 
we do not observe any convective motions, any grains 
spatial organization in the bulk and any grains order¬ 
ing close to the container wall. In addition, the driving 
based on freeze-thaw cycling is spatially homogenous. A 
structural analysis based on Vorono'i partition has been 
performed with numerical simulations. The distribution 
of the Vorono’i has shown that the packings are fully ran¬ 
dom during the whole process. Moreover, the packing 
fraction r) = V gra i n / /i v and the Vorono'i volume standard 
deviation a v are the main parameters characterizing the 
piles. 
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